package scu.maqiang.elasticity;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2Q42;
import scu.maqiang.mesh.Mesh2Q4;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.*;

public class PlateHole2DQ42 {

	public static void main(String[] args) {
	    double E = 3.0e5, Nu = 0.31e0;
	    String dir = "PlateHole2DQ\\";
	    String meshFile = dir + "PlateHole2DQ.msh";
	    Mesh2Q4 mesh2D = new Mesh2Q4();
		mesh2D.readMesh(MeshFormat.msh, meshFile);
	    int dispLabel = 20;
	    ScalarFunc dispRegion = (xy, label, param) -> xy[0] < 0.00001 ? 1.0 : -1.0;
	    mesh2D.setBoundaryLabel(dispRegion, null, dispLabel);
	    int forceLabel = 50;
	    ScalarFunc pressureRegion = (xy, label, param) -> xy[0] > 399.99999? 1.0: -1.0;
	    mesh2D.setBoundaryLabel(pressureRegion, null, forceLabel);
	    mesh2D.saveMesh(MeshFormat.msh, dir + "PlateHole2DQ.msh");
	    
	    FES2Q42 fs = new FES2Q42(mesh2D);
	    SRMatrix K = new SRMatrix(fs.GetNdof());
	    double[] ENu = {E, Nu};
	    fs.assembleStiff(ENu, BVPType.PLANESTRESS, K);
	    K.sort();
	    double[] x = new double[fs.GetNdof()];
	    double[] RHS = new double[fs.GetNdof()];
	    System.out.println(fs.GetNdof());
	    //fs.assemblePressure(new double[] {-1.0}, BVPType.COMMON, RHS, forceLabel);
	    fs.assembleFlux(new double[] {1.0, 0.0}, BVPType.COMMON, RHS, forceLabel);
	    System.out.println(MVO.sum(RHS));
	    fs.applyBC_MBN(K, RHS, Direct.All, 0.0, dispLabel);
		DirectSSolver solver = new DirectSSolver(K);
		x = solver.CholeskySolve(RHS, 1);
//	    NewIterSSolver solver = new NewIterSSolver(K);
//	    solver.PCGSSOR(RHS, x, 1.5, 1);
	    
	    double[][] disp = new double[2][fs.GetNdof()/2];
	    fs.extractComponent(x, disp);
	    System.out.println("Displacement:");
	    System.out.println("u min: " + MVO.min(disp[0]) + " max: " + MVO.max(disp[0]));
	    System.out.println("v min: " + MVO.min(disp[1]) + " max: " + MVO.max(disp[1]));
	    mesh2D.toTecplot(dir + "PlateHole2DRDisp.dat", disp);
	    mesh2D.toTecplotDeform(dir + "DeformHoleQuarter.dat", disp, 3.0e4);

//	    allocate(ex(fs%nDof/2), ey(fs%nDof/2), exy(fs%nDof/2))
//	    ex = 0.0D0
//	    ey = 0.0D0
//	    exy = 0.0D0
//	   
//	    call fs.computeStrain(u, v, ex, ey, exy)
//	    write(*, *) "Strain: "
//	    write(*, *) "ex min : ", minval(ex), "max: ", maxval(ex)
//	    write(*, *) "ey min : ", minval(ey), "max: ", maxval(ey)
//	    write(*, *) "exy min : ", minval(exy), "max: ", maxval(exy)
//
//	    allocate(tx(fs%nDof/2), ty(fs%nDof/2), tz(fs%nDof/2), txy(fs%nDof/2))
//	    allocate(eqStress(fs%nDof/2))
//	    tx = 0.0D0
//	    ty = 0.0D0
//	    tz = 0.0D0
//	    txy = 0.0D0
//	    eqStress = 0.0D0
//
//	    call fs.computeStress(ex, ey, exy, tz, E, Nu, tx, ty, txy, tz, PLANE_STRESS)
//	    write(*, *) "Stress: "
//	    write(*, *) "tx min : ", minval(tx), "max: ", maxval(tx)
//	    write(*, *) "ty min : ", minval(ty), "max: ", maxval(ty)
//	    write(*, *) "txy min : ", minval(txy), "max: ", maxval(txy)
//
//	    call fs.computeEqviStress(tx, ty, txy, tz, eqStress)
//	    write(*, *) "Von Mises Stress:"
//	    write(*, *) "min : ", minval(eqStress), "max: ", maxval(eqStress)
//
//	    outputFile = trim(dir)//'PlateHole2DRDisp.dat'
//	    call mesh2D.exportSolution(outputFile, u, v, eqStress)
//	    outputFile = trim(dir)//'PlateHole2DRStrain.dat'
//	    call mesh2D.exportSolution(outputFile, ex, ey, exy)
//	    outputFile = trim(dir)//'PlateHole2DRStress.dat'
//	    call mesh2D.exportSolution(outputFile, tx, ty, txy)
//
//	    outputFile = trim(dir)//'PlateHole2DRDeformed.dat'
//	    call mesh2D.exportMeshDeform(outputFile, u, v, 3.0D4)
	}
}
